library(naniar) # for missing data viz
library(tidyverse) # for data wrangling
library(skimr) # for statistical summaries
library(DT) # for nice tables
library(mice) # for missing data imputation
set.seed(123) # for reproducibility
The data is a subset of the 2009 survey from BRFSS (Behavioral Risk Factor Surveillance System), an ongoing data collection program designed to measure behavioral risk factors for the adult population (18 years of age or older) living in households.
More information: https://www.rdocumentation.org/packages/naniar/versions/0.5.1/topics/riskfactors
skim(riskfactors)
| Name | riskfactors |
| Number of rows | 245 |
| Number of columns | 34 |
| _______________________ | |
| Column type frequency: | |
| factor | 18 |
| numeric | 16 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| state | 0 | 1.00 | FALSE | 52 | 6: 12, 34: 12, 25: 11, 12: 10 |
| sex | 0 | 1.00 | FALSE | 2 | Fem: 153, Mal: 92 |
| marital | 1 | 1.00 | FALSE | 6 | Mar: 131, Wid: 41, Div: 39, Nev: 24 |
| pregnant | 215 | 0.12 | FALSE | 2 | No: 28, Yes: 2 |
| education | 1 | 1.00 | FALSE | 6 | 4: 85, 6: 75, 5: 59, 3: 14 |
| employment | 0 | 1.00 | FALSE | 7 | 1: 92, 7: 80, 2: 21, 8: 18 |
| income | 0 | 1.00 | FALSE | 10 | >75: 46, 35-: 38, 50-: 32, 25-: 28 |
| veteran | 3 | 0.99 | FALSE | 5 | 5: 198, 3: 35, 4: 6, 2: 2 |
| hispanic | 2 | 0.99 | FALSE | 2 | No: 221, Yes: 22 |
| health_general | 0 | 1.00 | FALSE | 6 | Goo: 82, Ver: 56, Exc: 45, Fai: 40 |
| health_cover | 0 | 1.00 | FALSE | 2 | Yes: 217, No: 28 |
| provide_care | 3 | 0.99 | FALSE | 2 | No: 181, Yes: 61 |
| activity_limited | 3 | 0.99 | FALSE | 2 | No: 170, Yes: 72 |
| drink_any | 2 | 0.99 | FALSE | 2 | No: 132, Yes: 111 |
| smoke_100 | 2 | 0.99 | FALSE | 2 | No: 126, Yes: 117 |
| smoke_days | 128 | 0.48 | FALSE | 3 | Not: 84, Eve: 26, Som: 7 |
| smoke_stop | 212 | 0.13 | FALSE | 2 | Yes: 19, No: 14 |
| smoke_last | 161 | 0.34 | FALSE | 6 | 7: 59, 5: 9, 6: 6, 8: 5 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| age | 0 | 1.00 | 58.11 | 17.50 | 7.00 | 48.00 | 59.00 | 70.00 | 97.00 | ▁▃▇▇▃ |
| weight_lbs | 10 | 0.96 | 174.27 | 44.70 | 96.00 | 144.00 | 170.00 | 195.00 | 410.00 | ▆▇▂▁▁ |
| height_inch | 2 | 0.99 | 66.35 | 3.80 | 57.00 | 64.00 | 66.00 | 69.00 | 75.00 | ▂▇▇▇▃ |
| bmi | 11 | 0.96 | 27.78 | 6.57 | 14.21 | 23.35 | 26.66 | 31.32 | 55.72 | ▃▇▃▁▁ |
| children | 0 | 1.00 | 0.42 | 0.88 | 0.00 | 0.00 | 0.00 | 0.00 | 4.00 | ▇▁▁▁▁ |
| health_physical | 0 | 1.00 | 4.12 | 9.11 | -9.00 | 0.00 | 0.00 | 3.00 | 30.00 | ▁▇▁▁▁ |
| health_mental | 0 | 1.00 | 3.14 | 8.19 | -9.00 | 0.00 | 0.00 | 1.00 | 30.00 | ▁▇▁▁▁ |
| health_poor | 113 | 0.54 | 5.46 | 10.76 | -9.00 | 0.00 | 0.00 | 5.00 | 30.00 | ▁▇▁▁▂ |
| drink_days | 134 | 0.45 | 9.29 | 10.20 | -9.00 | 2.00 | 4.00 | 15.00 | 30.00 | ▁▇▃▃▂ |
| drink_average | 135 | 0.45 | 1.41 | 2.57 | -9.00 | 1.00 | 1.00 | 2.00 | 9.00 | ▁▁▇▆▁ |
| diet_fruit | 8 | 0.97 | 324.71 | 338.11 | -9.00 | 60.00 | 209.00 | 365.00 | 2555.00 | ▇▁▁▁▁ |
| diet_salad | 8 | 0.97 | 175.48 | 140.00 | -9.00 | 52.00 | 156.00 | 261.00 | 730.00 | ▇▇▃▁▁ |
| diet_potato | 8 | 0.97 | 111.68 | 132.42 | -9.00 | 24.00 | 104.00 | 156.00 | 1095.00 | ▇▁▁▁▁ |
| diet_carrot | 8 | 0.97 | 91.52 | 101.32 | -9.00 | 24.00 | 52.00 | 156.00 | 365.00 | ▇▂▂▁▁ |
| diet_vegetable | 8 | 0.97 | 407.15 | 316.25 | -9.00 | 156.00 | 365.00 | 626.00 | 1825.00 | ▇▇▃▁▁ |
| diet_juice | 8 | 0.97 | 227.85 | 307.28 | -9.00 | 24.00 | 104.00 | 365.00 | 2555.00 | ▇▁▁▁▁ |
n_var_miss(riskfactors)
## [1] 24
gg_miss_upset(riskfactors)
The documentation on {naniar} says that the default option of gg_miss_upset is taken from UpSetR::upset - which is to use up to 5 sets and up to 40 interactions.
gg_miss_var(riskfactors)
gg_miss_var(riskfactors, show_pct = TRUE)
gg_miss_var(riskfactors,
show_pct = T,
facet = sex)
gg_miss_fct(x = riskfactors, fct = marital)
riskfactors %>%
group_by(marital) %>%
miss_var_summary() %>%
DT::datatable()
We want to assess the predictors of health_physical, a continuous variable corresponding to the number of days in the last 30 days for which the respondent said their physical health was good.
Explanatory variables of interest include the ones related to diet as well as health_poor. If we were to run a Complete Case Analysis regression analysis, we’d effectively be throwing away nearly half our sample.
riskfactors <- riskfactors %>% select(-c("pregnant", starts_with("smoke"), "drink_days", "drink_average"))
mice(). The default number of imputations is m = 5.quickpred() can be passed on to the pred argument if the researcher is interested in quickly selecting predictors in datasets that contain many variables. The threshold for correlation can be specified with mincor.Here we select predictors with a minimum correlation of \(\rho=.30\) :
impute <- mice(riskfactors,
m = 5,
pred = quickpred(riskfactors, mincor = .3),
print = FALSE)
Checking the class of impute
class(impute)
## [1] "mids"
plot(impute)
The plot shows the mean (left) and standard deviation (right) of the imputed values only. In general, we would like the streams to intermingle and be free of any trends at the later iterations.
As the documentation shows (https://www.rdocumentation.org/packages/mice/versions/3.5.0/topics/mice), many methods can be specified for imputation. To check which one our algorithm chose:
cbind(impute$meth)
## [,1]
## state ""
## sex ""
## age ""
## weight_lbs "pmm"
## height_inch "pmm"
## bmi "pmm"
## marital "polyreg"
## children ""
## education "polyreg"
## employment ""
## income ""
## veteran "polyreg"
## hispanic "logreg"
## health_general ""
## health_physical ""
## health_mental ""
## health_poor "pmm"
## health_cover ""
## provide_care "logreg"
## activity_limited "logreg"
## drink_any "logreg"
## diet_fruit "pmm"
## diet_salad "pmm"
## diet_potato "pmm"
## diet_carrot "pmm"
## diet_vegetable "pmm"
## diet_juice "pmm"
Looks like the algorithm used predictive mean matching (pmm) for the numerical variables, and either polytomous regression (polyreg) for unordered categorical variables or logistic regression for the binary categorical variables (logreg)
We can also change the method. Let’s use Bayesian linear regression - norm - for the bmi variable:
meth <- impute$meth
meth["bmi"] <- "norm"
cbind(meth)
## meth
## state ""
## sex ""
## age ""
## weight_lbs "pmm"
## height_inch "pmm"
## bmi "norm"
## marital "polyreg"
## children ""
## education "polyreg"
## employment ""
## income ""
## veteran "polyreg"
## hispanic "logreg"
## health_general ""
## health_physical ""
## health_mental ""
## health_poor "pmm"
## health_cover ""
## provide_care "logreg"
## activity_limited "logreg"
## drink_any "logreg"
## diet_fruit "pmm"
## diet_salad "pmm"
## diet_potato "pmm"
## diet_carrot "pmm"
## diet_vegetable "pmm"
## diet_juice "pmm"
We must re-run the imputation:
impute <- mice(riskfactors, meth = meth, print = FALSE)
We can plot the trace lines again to see what the convergence looks like:
plot(impute)
We’ve seen above that the lines do intermingle well, especially at the later iterations with just 5 iterations. Sometimes, it is useful to increase the number of iterations just to confirm that there is indeed no trend. We can increase the number of iterations to 20 by running 15 additional iterations using the mice.mids() function:
impute_20 <- mice.mids(impute, maxit = 15, print = FALSE)
plot(impute_20)
We can check our imputations for certain variables by comparing them to observed values - under the MCAR assumption, the imputations should indeed have the same distribution as the observed data. Under MAR, the distributions may be different but nevertheless very large differences warrant further investigation.
health_poor days:stripplot(impute, health_poor ~ .imp, pch = 20, cex = 2)
bmi:stripplot(impute, bmi ~ .imp, pch = 20, cex = 2)
bmi was imputed with Bayesian linear regression and (the range of) imputed values looks a little different than the observed values but there is still overlap.
Finally, we need to run the regression on each of the 5 datasets and pool the estimates together to get average regression coefficients and correct standard errors. The with() function in the mice package allows us to do this.
fit <- with(impute, lm(health_physical ~ age + bmi + sex + marital + children + employment + health_mental + activity_limited + provide_care + health_poor + health_general + diet_vegetable + diet_salad + diet_potato + diet_juice + diet_fruit + diet_carrot))
fit
## call :
## with.mids(data = impute, expr = lm(health_physical ~ age + bmi +
## sex + marital + children + employment + health_mental + activity_limited +
## provide_care + health_poor + health_general + diet_vegetable +
## diet_salad + diet_potato + diet_juice + diet_fruit + diet_carrot))
##
## call1 :
## mice(data = riskfactors, method = meth, printFlag = FALSE)
##
## nmis :
## state sex age weight_lbs
## 0 0 0 10
## height_inch bmi marital children
## 2 11 1 0
## education employment income veteran
## 1 0 0 3
## hispanic health_general health_physical health_mental
## 2 0 0 0
## health_poor health_cover provide_care activity_limited
## 113 0 3 3
## drink_any diet_fruit diet_salad diet_potato
## 2 8 8 8
## diet_carrot diet_vegetable diet_juice
## 8 8 8
##
## analyses :
## [[1]]
##
## Call:
## lm(formula = health_physical ~ age + bmi + sex + marital + children +
## employment + health_mental + activity_limited + provide_care +
## health_poor + health_general + diet_vegetable + diet_salad +
## diet_potato + diet_juice + diet_fruit + diet_carrot)
##
## Coefficients:
## (Intercept) age bmi
## 1.739e+00 -6.337e-02 1.503e-01
## sexFemale maritalDivorced maritalWidowed
## 2.048e+00 1.586e-01 1.593e+00
## maritalSeparated maritalNeverMarried maritalUnmarriedCouple
## -4.300e+00 -5.817e-01 -3.503e+00
## children employment2 employment3
## -6.484e-01 5.212e-01 -2.866e+00
## employment4 employment5 employment7
## -2.320e+00 1.163e+00 -4.259e-01
## employment8 health_mental activity_limitedNo
## -2.402e+00 1.047e-01 -4.212e+00
## provide_careNo health_poor health_generalVeryGood
## 4.186e-01 1.132e-01 4.179e-01
## health_generalGood health_generalFair health_generalPoor
## 6.289e-01 3.133e+00 1.679e+01
## health_generalRefused diet_vegetable diet_salad
## -8.066e+00 5.988e-04 1.030e-03
## diet_potato diet_juice diet_fruit
## -5.565e-03 4.084e-03 -4.019e-05
## diet_carrot
## -1.476e-03
##
##
## [[2]]
##
## Call:
## lm(formula = health_physical ~ age + bmi + sex + marital + children +
## employment + health_mental + activity_limited + provide_care +
## health_poor + health_general + diet_vegetable + diet_salad +
## diet_potato + diet_juice + diet_fruit + diet_carrot)
##
## Coefficients:
## (Intercept) age bmi
## 0.9990222 -0.0686284 0.1778563
## sexFemale maritalDivorced maritalWidowed
## 1.7289946 0.1417556 1.5580613
## maritalSeparated maritalNeverMarried maritalUnmarriedCouple
## -2.7448279 -0.4858848 -3.5094290
## children employment2 employment3
## -0.6290995 1.0762904 -2.5990223
## employment4 employment5 employment7
## -1.6601152 0.9011993 -0.8981470
## employment8 health_mental activity_limitedNo
## -2.9122112 0.1028946 -4.4034850
## provide_careNo health_poor health_generalVeryGood
## 0.2056059 0.1056196 1.1528270
## health_generalGood health_generalFair health_generalPoor
## 1.1309828 3.6589914 17.4273975
## health_generalRefused diet_vegetable diet_salad
## -6.9687402 -0.0000866 0.0043914
## diet_potato diet_juice diet_fruit
## -0.0041553 0.0038780 -0.0001935
## diet_carrot
## -0.0011747
##
##
## [[3]]
##
## Call:
## lm(formula = health_physical ~ age + bmi + sex + marital + children +
## employment + health_mental + activity_limited + provide_care +
## health_poor + health_general + diet_vegetable + diet_salad +
## diet_potato + diet_juice + diet_fruit + diet_carrot)
##
## Coefficients:
## (Intercept) age bmi
## 2.704129 -0.071143 0.135268
## sexFemale maritalDivorced maritalWidowed
## 1.830159 0.200566 1.503719
## maritalSeparated maritalNeverMarried maritalUnmarriedCouple
## -4.168948 -0.415450 -3.403838
## children employment2 employment3
## -0.749504 0.971755 -2.642386
## employment4 employment5 employment7
## -2.136584 0.841615 -0.552605
## employment8 health_mental activity_limitedNo
## -2.480937 0.104532 -4.453114
## provide_careNo health_poor health_generalVeryGood
## -0.030878 0.106614 0.647872
## health_generalGood health_generalFair health_generalPoor
## 0.988748 3.445597 16.724057
## health_generalRefused diet_vegetable diet_salad
## -7.735275 0.001092 0.002808
## diet_potato diet_juice diet_fruit
## -0.003753 0.004438 -0.001017
## diet_carrot
## -0.002032
##
##
## [[4]]
##
## Call:
## lm(formula = health_physical ~ age + bmi + sex + marital + children +
## employment + health_mental + activity_limited + provide_care +
## health_poor + health_general + diet_vegetable + diet_salad +
## diet_potato + diet_juice + diet_fruit + diet_carrot)
##
## Coefficients:
## (Intercept) age bmi
## 1.5677271 -0.0764285 0.1665056
## sexFemale maritalDivorced maritalWidowed
## 1.4635227 0.3138129 1.8212133
## maritalSeparated maritalNeverMarried maritalUnmarriedCouple
## -3.0856733 -0.7627401 -3.1658295
## children employment2 employment3
## -0.7966279 0.1802052 -2.6618967
## employment4 employment5 employment7
## -1.2421182 0.2962280 -1.0384293
## employment8 health_mental activity_limitedNo
## -2.8780833 0.0864878 -3.8797074
## provide_careNo health_poor health_generalVeryGood
## 0.3103232 0.1812182 1.3255393
## health_generalGood health_generalFair health_generalPoor
## 1.1020951 4.0242406 16.8904464
## health_generalRefused diet_vegetable diet_salad
## -5.9716591 0.0006650 0.0050499
## diet_potato diet_juice diet_fruit
## -0.0071418 0.0032890 -0.0006714
## diet_carrot
## 0.0009036
##
##
## [[5]]
##
## Call:
## lm(formula = health_physical ~ age + bmi + sex + marital + children +
## employment + health_mental + activity_limited + provide_care +
## health_poor + health_general + diet_vegetable + diet_salad +
## diet_potato + diet_juice + diet_fruit + diet_carrot)
##
## Coefficients:
## (Intercept) age bmi
## 2.6453979 -0.0764877 0.0868773
## sexFemale maritalDivorced maritalWidowed
## 1.8208563 0.3403289 2.0334635
## maritalSeparated maritalNeverMarried maritalUnmarriedCouple
## -6.8832296 -1.2065212 -4.8079140
## children employment2 employment3
## -0.6356130 0.7207182 -2.9229650
## employment4 employment5 employment7
## -2.0566007 0.3248994 -0.9867242
## employment8 health_mental activity_limitedNo
## -2.6273493 0.0906184 -3.4591495
## provide_careNo health_poor health_generalVeryGood
## 0.4729385 0.2205807 1.4180455
## health_generalGood health_generalFair health_generalPoor
## 2.0186421 4.4735674 16.6651833
## health_generalRefused diet_vegetable diet_salad
## -5.3054085 0.0007200 0.0037058
## diet_potato diet_juice diet_fruit
## -0.0071866 0.0039170 -0.0005845
## diet_carrot
## 0.0004970
The fit object contains the regression summaries for each data set. The new object fit is actually of class mira (multiply imputed repeated analyses). We can double-check:
class(fit)
## [1] "mira" "matrix"
This is the last step, where we pool the estimates from all 5 complete datasets to get average regression coefficients and correct standard errors.
pool.fit <- pool(fit)
summary_poolfit <- as.data.frame(summary(pool.fit))
summary_poolfit %>% DT::datatable()
class(pool.fit)
## [1] "mipo" "data.frame"